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The Big Raman spectrum of the two-dimensional S — 1/2 Heisenberg model is discussed within 
Loudon-Fleury theory at both zero and finite temperature. The exact T — spectrum for lattices 
with up to 6 x 6 sites is computed using Lanczos exact diagonalization. A quantum Monte Carlo 
(QMC) method is used to calculate the corresponding imaginary-time correlation function and its 
first two derivatives for lattices with up to 16 x 16 spins. The imaginary-time data is continued to 
real frequency using the maximum-entropy method, as well as a fit based on spinwave theory. The 
numerical results are compared with spinwave calculations for finite lattices. There is a surprisingly 
large change in the exact spectrum going from 4 x 4 to 6 x 6 sites. In the former case there is a single 
dominant two-magnon peak at u)/J » 3.0, whereas in the latter case there are two approximately 
equal-sized peaks at ui/J » 2.7 and 3.9. This is in good qualitative agreement with the spinwave 
calculations including two-magnon processes on the same lattices. The spinwave results for larger 
lattices show how additional peaks emerge with increasing lattice size, and eventually develop into 
the well known two-magnon profile peaked at uj/ J ~ 3.2 and with weight extending up to u)/ J ~ 4.6. 
Both the Lanczos and the QMC results indicate that the actual two-magnon profile is broader than 
the narrow peak obtained in spinwave theory, but the positions of the maxima agree to within a 
few percent. The higher-order contributions present in the numerical results are merged with the 
two-magnon profile and extend up to frequencies u/J « 7. The first three frequency cumulants of 
the spectrum are in excellent agreement with results previously obtained from a series expansion 
around the Ising limit. Typical experimental B\ g spectra for La2Cu04 are only slightly broader 
than what we obtain here. The exchange constant extracted from the peak position is J » 1400K, 
in good agreement with values obtained from neutron scattering and NMR experiments. We discuss 
the implications of our present results for more sophisticated theories of Raman scattering suggested 
recently. 

PACS numbers: 75.40.Gb, 75.40.Mg, 75.10.Jm, 75.50.Ee 



I. INTRODUCTION 



The magnetic properties of the parent compounds of 
the high-T c cuprate superconductors can be well ac- 
counted for by weakly coupJed two-dimensional (2D) 
Heisenberg antiferromagnetsll Neglecting the weak inter- 
layer coupling, the model is defined by the Hamiltonian 



/T = J^S 4 -S„ (J>0), 



(1) 



where Si is a spin-1/2 operator at site i on a square lat- 
tice and denotes a pair of nearest-neighbor sites. 
The most well studied amo«g|ih.e_antiferromagnetic lay- 
ered cuprates is La2Cu04,BHtai3'El with a Neel ordering 
temperature TV ~ 300 K. For T > TV, the temperature 
dependence of the spin-|Correlation length, measured us- 
ing neutron scattering,0 is in good agreement with that 



of a single-layer Heisenberg model with J 1500 K.ii 
The spinwave spectrum of the Heisenbergjpaodel is well 
reproduced over the entire Brilloin zone.erel The NMR 
relaxation rates l/Tj and I/T^g-, which probe the low- 
frequency spin dynamics n |also show remarkahlfi agree- 
ment between experimenttlcj and theory.BESOO 

In contrast to thes£_jsuccess stories, the experimen- 
tal Raman spectrumEjO'Ej shows significant deviatioBa 
from calculations for the 2D Heisenberg model.E3ll2Hl3't3 
Within this description of the CuC>2 layers, the stan- 
dard theory of Raman_scattering is based on the Loudon- 
Fleury (LF) couplingo between the light and the spin 
system. The coupling is obtained in second order pertur- 
bation theory with virtual states .containing one doubly 
occupied site, and is given byEJ'EII 



H 



LF 



(id) 



(Ej n ■ <7ij)(E out • <Tij)Si ■ Sj 



(2) 



1 



Here Ei n and E out are the polarization vectors of the 
incoming and scattered light and is the unit vector 
connecting sites i and j. In terms of the eigenstates {|n}} 
of the Heisenberg model, the frequency dependence of the 
scattering intensity at inverse temperature [3 is given by 
Fermi's Golden Rule: 

m 

J2\(n\H LF \m)\ 2 6(LU-[E n -E m }). (3) 

n 

Most theoretical work has focused on the scattering in the 
Big symmetry channel. This corresponds to E in along a 
diagonal of the square lattice, and E ou t perpendicular to 
Ei„. The Big coupling can thus be written as 

/// '• E s > - s - E s <- s - ( 4 ) 

where (i,j) x and (i,j) y denote links in the x and y di- 
rections, respectively. 

For La 2 Cu04, the Bi g spectrum has a broad asymmqt-. 
ric peak at u> ~ 3 J with a tail extending tows; 7 — 8J.E3 
In some cases there is a shoulder-like structure atwx 4 J. 
Within spinwave theory the Bi g scattering is dominated 
by two-magnon excitations. The two-magnon profile is 
peaked around u .Ki. 3J in good agreement with the 
experiments.E3'U'HL3 However, the large width of the 
experimental spectrum has not been reproduced within 
spinwave theory. The first three frequency cumulants of 
the spectrum have also been calculated using a series ex- 
pansion around the Ising limit ,E3 and are in reasonable 
agreement with the experimental values. The moments 
obtained in spinwave theory to order 1/S (two magnon 
excitations only) are in poor agreement with the series 
results which in principle include multi-magnon contri- 
butions to all orders. Unfortunately, the full frequency 
dependence is not accessible with the series expansion 
method. Exact diagonalization has been used tq-com- 
pute the exact LF Raman profile for small lattices.li3 For 
the 4x4 system there is a single dominant two-magnon 
peak at uo/J = 2.98. Weight present at u> « 5 J has 
been attributed to four-magnon processes, but its relative 
strength is much smaller than the weight found experi- 
mentally in this frequency region. The tail at higher fre- 
quencies is absent. Despite this, the first three frequency 
cumulants are in approximate agreement with both ex- 
periments and the secies expansion results. 

Canali and GirvirJIZI carried out a spin-wave expansion 
including also four-magnon excitations (which enter in 
order S~ 2 ). The narrow width of the two-magnon peak 
was found to be stable with respect to inclusion of the 
higher-order processes. The relative contribution from 
four-magnon states in this calculation is less than 3%. 
The high energy of the four-magnon weight nevertheless 
leads to first and second cumulants that are much closer 
to those obtained in the series expansion and exact di- 
agonalization studies. The spinwave result for the third 



cumulant is, however, significantly larger than the series 
expansion value. It was argued that this is due to inter- 
actions between four magnons that were neglected in the 
1 / S 2 calculation, and that the relative four-magnon con- 
tribution must h& w 10% in order to reproduce the first 
three moments £3 The conclusion that there is less high- 
energy weight than in typical experimental spectra then 
still remains. However, the apparent inability of a very 
sophisticated spinwave calculation to fully capture the 
four-magnon processes raises some concerns about this 
approach for calculating the actual line shape. Further- 
more, Chubukov and Frenkel have recently questioned 
the stage at which the spin S was set to 1/2 in the pre- 
vious spinwave calculations. They kept S large and car- 
ried out an expansion of the profile around its peak posi- 
tion before expanding in 1/ S and evaluating the result at 
S = 1/2. The two-magnon profile obtained this way has 
a width almost three times larger than the "standard" 
one, and the second frequency cumulant is therefore in 
better agreement with the series result. However, the 
agreement with the first cumulant is actually worse, due 
to the significantly larger low-frequency weight. Hence, 
this result also has to be viewed with some caution. 

Experimentally, sLg-Hificant scattering is also observed 
in the Ai g channel. t2l With the standard LF coupling 
inelastic scattering in this symmetry is not possible for 
the Heisenberg model with only nearest-neighbor inter- 
actions, since the x and y terms in Eq. (^) are added 
in this case and Hpp then commutes with the Hamil- 
tonian. Adding a next-nearest-neighbor term to the LF 
coupling leads, to Ai g scattering, but has no effect in the 
Big channel.E3 Based on the frequency moments obtained 
in the series expansion, this was argued to be a mech- 
anism that could explain both the Bi g and Ai g spec- 
tra. However, the narrow Bi g line-shape obtained in 
other calculations remains an unresolved issue for this 
scenario. Including a next-nearest-neighbor term ~ J2 
in the Hamiltonian would also lead to Ai g scattering, 
and probably a broadened Bi g spectrum. The relatively 
large Ai g intensity seen experimentally would then likely 
require a larger J2 than allowed by other experiments, 
but detailed calculations have not been carried out within 
this model. Other interactions, such as the so called four- 
spin cyclical exchange ,E£rt3 have also been suggested to 
account for the differences between theory and experi- 
ment. Analytical calculations as detailed as those for the 
standard LF theory (only nearest-neighbor interactions 
in both the Hamiltonian and the Raman operator) have 
also not yet been carried out within these theories. 

It is clear that LF theory is not sufficient to capture 
all aspects of Raman scattering in layered cuprates. For 
example, resonant scattering (occurring when the fre- 
quency of the incoming light is comparable to the charge- 
transfer gaplxan of course not occur within the Heisen- 
berg model.Efl However, focusing on nonresonant Bi g 
scattering, it is not yet clear what the actual line shape 
is within LF theory. As discussed above, the two key 
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questions of the width of the dominant two-magnon pro- 
file and the relative weight of the higher-order contribu- 
tions above the two-magnon cut-off remain incompletely 
answered within spinwave theory. The exact diagonal- 
ization studies carried out so far are limited to lattices 
too small for reliable quantitative extrapolation to the 
thermodynamic limit. There is hence a definite need 
for accurate non-perturbative numerical calculations for 
larger lattices. Conclusive results would provide a more 
solid basis for estimating effects not included in the LF- 
Heisenberg model. It should also be noted that LF theory 
is the standard framework in which Raman scattering has 
been interpreted also in several othep4ow-dimcnsional an- 
tiferromagnetic 5=1/2 systems.E3'E3 A satisfactory res- 
olution of the 2D Heisenberg case would therefore be of 
more general interest as well. Finally, knowing the exact 
Big Raman spectrum should help to shed light on the 
applicability of spinwave theory for calculations of dy- 
namic processes involving excitations of more than one 
magnon. 

Here we report new exact diagonalization results for 
systems with up to 6 x 6 spins, which is the largest cur- 
rently accessible with this method. We also obtain ap- 
proximate spectra for up to 16 x 16 spins by Maximum- 
Entropy (Max-Ent) analytic continuation of quantum 
Monte Carlo (QMC) data^rUsAng the stochastic series 
expansion QMC techniqucOEH we have calculated the 
imaginary-time LF correlation function and its first two 
derivatives at temperatures low enough to give ground 
state results. The derivatives are used as supplementary 
information in the analytic continuation. We also con- 
sider a more phenomenological approach of fitting the 
two-magnon profile obtained within spinwave theory to 
the imaginary-time data, adding a Gaussian at higher fre- 
quency to model the higher-order contributions. In order 
to study the effects of temperature, we apply the QMC 
+ Max-Ent methods also at non-zero temperatures. 

We find that the B\ g Raman spectrum of the 6x6 
lattice has two dominant peaks at oj / J ~ 2.7 and 3.9, in 
sharp contrast to the single dominant peak at to/ J = 2.98 
previously found for the 4x4 lattice. These results are 
qualitatively reproduced within spin wave theory includ- 
ing only two-magnon excitations and magnon-magnon in- 
teractions (treated within an RPA scheme). Spinwave 
results resembling the infinite-size two-magnon profile 
are seen only for much larger lattices. Using the QMC 
data, the first three frequency cumulants in the ther- 
modynamic limit can be reliably estimated. They are 
in excellent ..agreement with the previous series expan- 
sion results. E3 Both the Lanczos and the QMC results 
indicate that the dominant peak at u> ~ 3.3 J is slightly 
broader than the standard two-magnon profile, but not 
as broad as the the one recently obtained by Chubukov 
and Frenkelo We find no evidence for a gap between 
the two-magnon profile and the higher-order contribu- 
tions. They appear to be completely merged together 
and extend up to uj ~ 7 J. 

Finite-temperature results for T/J < 0.25 are very 



similar to the ground state results. For higher tempera- 
tures there is a significant growth of the low-frequency 
spectral weight, as also found in a previous exact di- 
agonalization study of a 4 x 4 lattice.EJ We find that 
the temperature at which this effect becomes significant 
decreases with increasing system size, due to the large 
finite-size gaps present in the smaller systems. 

Over-all, our results are in closer agreement with ex- 
periments than previous exact diagonalization and spin- 
wave calculations, but we conclude that the Heisenberg- 
LF spectrum is nevertheless not quite as broad as typical 
experimental spectra. We argue that our results give 
more credibility to. pypraosed broadening mechanisms in- 
volving phononsJlil'EirESH'Eil 

The rest of the paper is organized as follows. In Sec. II 
we review various spinwave calculations of the B\ g Ra- 
man profile. We also present results for small lattices, 
which are compared with exact diagonalization spectra 
in Sec III. Our T — results from QMC simulation and 
numerical analytic continuation are discussed in Sec. IV. 
Effects of finite temperature are considered in Sec. V. In 
Sec. VI we summarize and discuss our results and impli- 
cations for mechanisms proposed to lead to a broadening 
of the B\g profile. In an Appendix we present some tech- 
nical details of the QMC calculation of the imaginary- 
time correlation function and its derivatives. 



II. SPINWAVE THEORY 

The Raman spectrum can be easily computed in the 
spinwave approximation. Various improvements of the 
linear spinwave calculation can be also applied to the 
Raman scattering amplitude. For example, the resid- 
ual interactions between spin waves, which play a cru- 
cial role in the Raman excited states, can be included 
at the RPA level. Here we discuss linear spin wave the- 
ory and the effect of magnon-magnon interactions on the 
B\g spectrum. The primary purpose of the calculations 
discussed here is to qualitatively understand the effects 
of finite size, which will be important for interpreting 
the numerical results presented in the following sections. 
More sophisticated spinwave calculations including the 
quantum fluctuations of the ground state, as well as fi- 
nal states with four magnons, have been carried out be- 
fore, as discussed in the Introduction. Here we note some 
problems with the analytical calculations of the Raman 
profile which motivate our renewed efforts to obtain ac- 
curate non-perturbative numerical results. 

In the antiferromagnetic ground state we are consider- 
ing, a bosonic representation of the spin operators can be 
introduced on each sub latt ice by using the usual Dyson- 
Maleev transformational^. On sublattice A, the trans- 
formation reads 

S* = S-a\a t (5a) 
S,+ = V25(l-^)a,: (5b) 
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s; 



2Sa 



(5c) 



where a[ creates a boson, i.e., a spin-1 magnon, at site i. 
Similarly, on sublattice B, 



5T = 



-s 



I— t / b % 



2Sbi 



(6a) 

(6b) 
(6c) 



As usual, as long as small fluctuations around the Neel 
ordered phase are considered, one keeps only quadratic 
terms in the Hamiltonian (which are the dominant terms 
in a 1/5 expansion). Therefore, in this approxima- 
tion, the Hamiltonian can be diagonalized by a Bogoli- 
ubov transformation in terms of spinwave excitations (or 
magnons); 



/?k = "k&k 



where the coefficients (> 0) are given by, 



1 ( 1 



1 f 1 



(7a) 
(7b) 



(7c) 
(7d) 



We have defined 7k = J2s e /Z as a sum over the Z 
nearest neighbours of the site at the origin. In our case 
(square lattice), 7k = (cosfc x + cosfc a )/2. The resulting 
well known linear spinwave Hamiltonian reads, 



Hsw = E Q +Y] Wk(a k a k + /^AJ, 



k 



(8) 



with the dispersion relation Wk = JSZy/l — 7^. Due 
to the decomposition into two sublattices, the reciprocal 
space is reduced to the magnetic brillouin zone (MBZ). 

The LF operator can be easily expressed as a quadratic 
form in terms of the spin wave operators. At zero tem- 
perature the ground state contains no bosons if magnon- 
magnon interactions are neglected so that, in a first ap- 
proximation, one only keeps constant terms or terms in- 
volving the creation of magnons, giving 



H 



LF 



-aNS 2 E out ■ E in (9) 

2aS (EontK [cos k x (ul + t£) - 2u k «k] (10) 

k 

E v out El[co S k v (ul + vl) - 2u k v k ]) 4/3 k . (11) 



(12) 



In the Big geometry, the matrix element is thus 



(f\H R \i) 



cos k x — cos k y 



The Raman intensity obtained from Fermi's Golden Rule, 
Eq. (|), is then 



J( W ) ex 



{c0S \Z™{ ky)2 S("-™*), (13) 



where ilk = ^JSy/l — 7^ is the frequency of the magnon. 
This expression exhibits a divergence at u> = 8J5 since 
the density of states diverges at the boundary of the 
MBZ. 

It is well known that this result is strongly modified 
when one takes into accou nt , the .. rn agnon- magnon in- 
teractions in the final stateOOEZlEJ In order 1/5, the 
Dyson-Maleev transformation generates in the Hamil- 
tonian quartic terms in the bosons operators. One 
way of treatine^this interaction is to keep the term 
a k<^k^-p a p'— ' w hi crl is responsible for multiple scat- 
tering of two magnons in the vacuum. This part gives the 
dominant contribution to the magnon-magnon scattering 
coming from the region near the MBZ boundary where 
the density of states diverges. Further simplification re- 
sults from the vanishing of 7k at the MBZ boundary so 
that it is reasonable to replace 7k by zero fcp^all k. This 
leads to an effective interaction of the form,E-3 



H ^ = ~J[zZzZ fk-poi^-p ^- 

k P 



(14) 



Following Refs. ||, and expressing 7k- P as a func- 



tion of its symmetric terms 7k-p = 7k + 7 P + + 7k 7p 
7k + 7 P + +7r7p~ 



7k = (cos k x ± cos ky) /2 
7£ ± = (sin k x ± sin k v )/2, 



(15) 
(16) 



it can be shown that multiple diffusion RPA series only 
contains terms involving 7p~ factors. The final RPA ex- 
pression for the Raman intensity is given bypJ 



I(to) oc Im 



R(u) 



l + i?(w)/45 



with 



8J5 



E 



(cos k x 



cos k y ) 2 



2^L 



■ is 



(17) 



(18) 



In the thermodynamic limit, Eq. ( |l8| ) for 5 = 1/2 leads 
to a narrow two-magnon peak around u> — 2.78J which 
extends up to to/ J = 4. In order to be able to directly 
compare spinwave results with exact spectra obtained 
with the Lanczos diagonalization method and approxi- 
mate results of QMC and Max-Ent analytic continuation 
(presented in the following two sections), we have also 
evaluated (|l7|) for small lattices. Results for L x L clus- 
ters with L = 4, 6, 8 and 10 are shown in Fig.Jj], along 
with the corresponding non-interacting form (|l3|). It is 
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clear that the continuous absorption band for infinite size 
is obtained from the accumulation of a series of peaks for 
increasing cluster sizes. However, the number of peaks is 
still very small even for a lattice of 100 sites. 

In the above calculation the magnon-magnon interac- 
tions have only been included in the final two-magnon 
state. The main effect of interactions in the ground 
state is to renormalize the spinwave velocity; c — > Z c c, 
where to order 1/S the renormalization factor Z c = 1.158 
and arises solely feopi the normal ordering of the quar- 
tic magnon terms £3 Hence, in a phenomenological way, 
the 1/5* corrections can be taken into account by shifting 
the energy scale by a factor Z c , leading to a B\ g profile 
peaked at co/J = 3.22. However, Canali and Girvin have 
shown that the renormalization in fact enters in a non- 
trivial way in the expression for the -Raman profile, and 
hence there are other effects as welLED Nevertheless, the 
end result for S — 1/2 does not differ much from Eq. (|T 
with Z c = 1.158, as shown in Fig. [|. 

Chubukov and Frenkel have recently raised questions 
about the stage at which one should set S = 1/2 in the 
spinwave calculation. They argued that one should first 
expand the large-S expression for the profile ( |l7| ) around 
its peak position, and only after that set S = 1/2.CJ The 
position of the maximum then remains approximately 
the same, but the profile is considerably broader, as also 
shown in Fig. ^. The better agreement r-with the fre- 
quency moments calculated by Singh et al23 was used as 
support in favor of the broader peak. However, it should 
be noted that the agreement with the first cumulant is 
actually significantly worse, due to the much slower de- 
cay of the weight on the low-frequency side of the peak. 
In fact, as seen in Fig. ^, the Chubukov- Frenkel profile 
extended towards lower frequencies in a way similar to 
the experimental spectrum, but the high-frequency tail 
is still of course missing. The high-frequency scattering 
was argued to be dominated by resonant scattering in 
typical experimental situations, and would hence not be 
explainable by LF theoryO 

To higher orders in 1/S, states enter in which 4, 6, 
e.t.c, magnons are excited (the true ground state is a 
linear combination of states containing any even number 
of magnons, and the Raman operator can create or anni- 
hilate one or two pairs of magnons). Canali and Girvin 
included four-magnon excitations (order 1/S 2 ) but np=. 
glected interactions involving more than two magnons .til 
The two-magnon profile obtained this way is very simi- 
lar to the Canali-Girvin 1/S result discussed above. The 
relative contribution from four-magnon processes is less 
than 3%, bu±-,is likely strongly affected by the neglected 
interactions .O In fact, although the small four-magnon 
contribution is sufficient (because of its rather high en- 
ergy) to change the first and second frequency cumulants 
to values in close agreement with the series expansion 
results, the third cumulant remains far off. It was ar- 
gued that this inconsistency is due to the neglected in- 
teractions among four magnons (which may even lead to 
bound states), and that such interactions would bring 



the four-magnon peak position down in frequency.0 The 
relative four-magnon weight would then have to increase 
to w 10% in order to satisfy the first three moments. In 
our opinion, this rough estimate indicates that the ap- 
proximations made in the 1/S 2 calculation may in fact 
be serious. In particular, it is not clear that the two- 
magnon and four-magnon contributions will be well sep- 
arated from each other if the four-magnon weight moves 
down and increases by a factor of 3 or more. This, in 
turn, may lead to considerably stronger interference ef- 
fects that may cause changes also to the upper edge of 
the two-magnon profile (which then no longer would arise 
from two-magnon excitations only). 

There are hence two major concerns with the spinwave 
calculations that have to be addressed: 1) The stage at 
which S is set to 1/2, leading to two very different two- 
magnon line shapes. 2) The contributions from processes 
including four or more magnons, which are very difficult 
to capture completely within spinwave theory. In this 
situation it is clearly useful to consider non-perturbative 
numerical methods. We discuss two complementary ap- 
proaches in the next two sections. 



III. EXACT DIAGONALIZATION 

In this section, we compute the exact Raman spectrum 
on clusters with up to ./V = 36 sites by use of the Lanczos 
diagonalization algorithm. In this approach, the Raman 
spectrum is obtained from a continued fraction, 

J(u) = --Tm({Q\Hl F — |- ^H LF \0)\, (19) 

where |0) is the ground state of energy Eq which can be 
easily calculated with the Lanczos method, and e is a 
small imaginary part added to give a finite damping of 
the ^-functions. 

Results for several square and tilted lattices are shown 
in Fig. H| along with the RPA spinwave results discussed 
in the previous section. We have shifted the spinwave 
results by the renormalization factor Z c ~LJ8 obtained 
using several different numerical methodsEj'O (this value 
is also in clpae agreement with the 1/S 2 spinwave value 
Z c = 1.177113). Spin wave theory clearly correctly pre- 
dicts the number of the dominant peaks, which hence can 
be characterized as two-magnon peaks. There are, how- 
ever, some discrepancies in the peak positions and their 
relative weights. Most notably, for the largest lattice 
(6x6), the separation between the two peaks in the ex- 
act spectrum exceeds by a factor of more than 1.5 that of 
the spinwave result. This may well be an indication that 
the correct two-magnon profile is broader than, the stan- 
dard profile obtained with spinwave theory.EEnlZl Whether 
or not it is as broad as that obtained by Chubukov and 
Frenkel (see Fig. ||) cannot be determined from the re- 
sults for these small lattices, however. We will return to 
this important issue in the next section. 
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It can be noted that for the 4x4 lattice there is a small 
peak at u / J ~ 4.5 both in the spinwave and the Lanczos 
results. This then suggests that it is a two-magnon peak, 
in contrast to previous claims that it arises from a four- 



magnon state.E3 For the larger lattices there is visible 
weight extending up to uj/J w 7, which is not present 
in the spinwave results and hence is due to processes 
involving more than two magnons. For the 6x6 lattice 
the relative weight of these contributions is about 10%. 

It is difficult to scale the full shape of the spectrum to 
infinite system size. The first few frequency cumulants 
can be expected to converge rather quickly, however, and 
have also previously been calculated using a series ex- 
pansion method as already discussed. The nth frequency 
moment is given by (at T = 0) 



p n = I dujuj n I (u>) . 
Jo 

The first cumulant Mi = pi, and for n > 1, 

(M n ) n = - f°° du(u ~ Pi) n I(u). 
Po Jo 



(20) 



(21) 



The results for 6 x 6 (4 x 4) are Mi = 3.524 (3.244), 
M 2 = 0.8686 (0.797), and M 3 p 0.9576 (1.141). The 
previous series expansion results^ 2 ] are Mi = 3.58 ± 0.06, 
M 2 = 0.81 ± 0.05, and M 3 = 1.00 ± 0.14. Hence, the 
6x6 cumulants show an improved and good agreement 
with the series results. However, since there are signifi- 
cant differences between the 4x4 and the 6x6 lattices, 
the results may still change in the thermodynamic limit. 
Unfortunately, using also the results for TV = 26 and 32, 
the data do not fall on smooth curves (see Fig. ^ in the 
next Section) , and it is not possible to extrapolate to the 
thermodynamic limit using only these Lanczos results. 
In the next Section we will calculate the cumulants for 
much larger systems using QMC data. 



IV. QUANTUM MONTE CARLO AND 
MAXIMUM-ENTROPY ANALYTIC 
CONTINUATION 

Real-frequency dynamic properties cannot be obtained 
directly using QMC methods. Instead, the correspond- 
ing imaginary-time dependent correlation function has 
to be calculated, and numerically continued to real fre- 
quency. For the Raman spectrum defined by Eq. (||), the 
imaginary-time function is given by 



G(r) = (Hlf(t)Hlf(0)), 



(22) 



where Hlf{t) = e T H^pe 



-tH 



The analytic continua- 



tion to real frequency amounts to inverting the integral 
relation 



G(r) 



dcuI(Lo)e 



(23) 



With G(t) obtained only to within a statistical error 
from a QMC simulation, the spectrum I(cu) cannot be 
uniquely determined— In the Max-Ent approach to this 
difficult problem,oEJ a unique solutions is defined as 
that minimizing 



Q = ]pC 2 - aS, 
where S is the entropy of the spectrum, 

S -I dujl(uj) In [I(w)/m(w)], 



(24) 



(25) 



defined with respect to a "default" model m (both / and 
m are here assumed to be normalized to unity). G(r) 
is calculated for a discrete set of times Tj. A given 
I(cj) corresponds to unique values of G(Tj) according 
to Eq. (^3|). The deviation from the actual calculated 
GQMc(i~i) is quantified by x 2 . Since the statistical er- 
rors <7j of GQMc{ T i) at different times are correlated (see 
Fig. [l3| in Appendix A), \ 2 should be defined in terms of 
the inverse of the full covariance matrix G, 

X 2 = J2[G(n) - Gqmc ^ij 1 [G( T j) - G Q mc{t 3 )} • 

(26) 

We here parametrize the spectrum in terms of ~ 
200 — 400 equally spaced delta-functions 6(ui — uii) for 
lu, > 0: 



E 



Ii8{tjJi - lo), 



(27) 



and a smooth continuous spectrum is then represented 
by the curve connecting the amplitudes U (or by giv- 
ing the 5-functions a width of the order of the frequency 
spacing, which gives a very similar curve) . The negative 
part of the bosonic spectrum is given by detailed bal- 
ance: I{—oS) = e~P u I(u>). We use a flat default model 
for oj > 0. The parameter a in Eq. ( pi| ) is determined 
iteratively so as to satisfy the "classic" Max-Ent crite- 
rion, resulting in (within the assumptions of the Max-Ent 
method) the spectrum with the highest probability given 
the QMC data. 

For calculating C?(tJ gR use the stochastic series ex- 
pansion QMC method,EZrE3 as discussed in Appendix A. 
With this method, derivatives of G(r) can also be di- 
rectly calculated. We here use the first two derivatives 
as supplementary information in the Max-Ent method. 
The nth derivative of G(t) is related to I(u>) according 
to 



(-1)" 



dujuj n l(uj)e 



(28) 



It is a straight-forward matter to modify the Max-Ent 
procedures to include also the first few (in our case 2) 
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of these derivative relations in addition to the original 
analytic continuation equation (^3|). The use of deriva- 
tives was first suggested by Schuttler and Scalapino in 
their pioneering work on numericaL-ajialytic continuation 
based on \ 2 fitting to QMC dataJl3 To our knowledge, 
the Max-Ent method has not previously been used with 
derivative information. It should be noted that the nth 
frequency moment p n is given by the r = derivative: 



p n = (-iy 



t GW(r^0) 
G(r^O) ' 



(29) 



Enforcing known frequency moments has been previously 
used to improve the resolution of the Max-Ent method.CJ 
The derivative information goes beyond this by enforcing 
also the "moments" defined with r > in Eq. (p9|). The 
derivatives can of course be expected to improve on the 
Max-Ent procedure only if they can be calculated accu- 
rately enough to contain information not already present 
in the calculated G(r). Typically, the statistical errors 
increase with increasing derivative order n. In our case, 
the first two derivatives appear to be useful, although 
spectra obtained with only G(t) are not dramatically 
different. 

Next, we present results for systems of size Lx L, with 
L = 4, 6, 8, and 10. Although considerably larger lat- 
tices can be studied with the QMC method, the physical 
information we are interested in here requires very ac- 
curate results for G(t). It is therefore more appropriate 
to concentrate the computational resources on obtaining 
reliable results for moderate system sizes. Comparing re- 
sults for L = 4 — 10 should also be sufficient for making 
statements about the thermodynamic limit. We also car- 
ried out some simulations for L = 16, but the statistical 
error are significantly larger in this case and the continua- 
tion to real frequency is therefore less reliable. In order to 
obtain ground state results, the simulations were carried 
out at inverse temperatures as high as j3 = 8L. Results 
obtained with (3 = AL are indistinguishable within sta- 
tistical errors, indicating that contributions from excited 
states indeed are negligible at these low temperatures. 

We begin by showing in Fig. || our results for the 
logarithm of the normalized imaginary-time correlator 
gir) = G(t)/G(0). In the same figure we also show the 
relative statistical error, er rc i(r), of g(r). Since the re- 
sults for all the system sizes have comparable errors, one 
can expect the Max-Ent continuation to real frequency 
to resolve structure on roughly the same scale. Already 
from this imaginary-time data it is clear that the real- 
frequency spectum has dominant weight at u> f=a 3J for 
all system sizes, as In [g(r)] decays approximately linearly 
with r in a sizable regime, with slope « —3. For the 
larger systems a slight upward curvature can be noted, 
indicating that there is spectral weight also below 3 J. 

We find that the shape of the Raman spectrum ob- 
tained with the Max-Ent method is very sensitive to 
the statistical fluctuations in the QMC data. Carrying 
out the Max-Ent procedures with different subsets of the 



available imaginary-time data always gives a dominant 
peak close to uj/J = 3, but the peak width and asym- 
metry show large variations. We therefore consider it 
appropriate to define the spectrum corresponding to the 
full set of imaginary time data as an average over suit- 
ably defined subsets. For this purpose we use the so 
called bootstrap methodo in the following way. 

With the simulation data for some quantity A divided 
into M "bin averages" Ai in the standard way, a boot- 
strap sample Ab is defined as 



A, 



1 M 



(30) 



whre Ri is a randomly chosen bin (i.e., the number of 
bins chosen is the same as the total number of bins, allow- 
ing, of course, for multiple selections the bins). Since the 
Max-Ent procedure is highly non-linear, the average over 
a large number of separately Max-Ent continued boot- 
strap samples of imaginary-time correlation functions can 
be different from the continued full average. We argue 
that the bootstrap average is more meaningful since sta- 
tistical fluctuations are averaged out considerably. 

In Fig. |5| we show Max-Ent results for 10 bootstrap 
samples of 4 x 4 QMC data. All the spectra have a dom- 
inant peak very close to the correct position oj/J = 2.98, 
as well as a structure at higher frequency. There are, 
however, very large variations in the peak width and in 
the position of the high-frequency weight. The average 
over 500 bootstrap samples is shown in Fig. ^. The exact 
Lanczos result with a damping e/J = 0.1 is quite well 
reproduced, except that the small peak at w/J = 4.5 
is not present. It can be noted that the main peak of 
the average spectrum is narrower than most of the "typ- 
ical" bootstrap samples (see Fig. ||), contrary to what 
might have been expected. This is clearly due to the fact 
that the position of the peak shows very small variations 
compared to the variations in the peak width and that 
for some bootstrap samples the peak is very sharp. 

Results for the larger lattices are also shown in Fig. ||. 
In the 6x6 spectrum the two main peaks are clearly 
resolved. The weight present at higher frequency cannot 
be resolved as a separate structure, however, and instead 
causes the shift by « 15% of the second peak. As the sys- 
tem size grows, the number and density of peaks increase, 
and only a single structure can then be resolved. For the 
8x8 lattice, the spinwave result shown in Fig. [j] has only 
one dominant peak. The Max-Ent result for this size is, 
however, very broad, indicating that the relative weight 
distribution among the peaks obtained in spinwave the- 
ory is not reliable (signs of this is seen also in the exact 
6x6 spectrum in Fig. ||). In particular, the Max-Ent 
spectrum has much more low-frequency weight. This is 
the case also for the 10 x 10 lattice. The spectrum has 
a more pronounced peak than for 8x8, indicating that 
the individual ^-functions begin to group into a profile 
peaked around lo/J~ 3.5. It should be noted that the 
procedures we are using can be expected to work better 



7 



for the larger systems, for which the distribution of 6- 
functions are better approximated by a single continuous 
structure. 

In Fig. fj] we show the results for the short-time behav- 
ior of the ratio G' b '(t)/G(t), along with the correspond- 
ing curves obtained from the Max-Ent results. According 
to Eq. (p9|), the first two frequency moments can be di- 
rectly obtained from the r = points. The first moment 
can be accurately extracted this way. In the case of the 
second moment the statistical fluctuations grow large as 
t — > 0, but the extrapolation provided by the Max-Ent fit 
still gives a quite stable result. We also extract the third 
moment from the Max-Ent spectra. For both L — 4 and 
L = 6 the results are in excellent agreement with the ex- 
act results obtained with the Lanczos method in Sec. III. 

Fig. U shows the system size dependence of both the 
QMC and the Laaczos results for the cumulants, along 
with the previousE3 infinite-size series results by Singh et 
at. We also include the first and second cumulants ob- 
tained for a 16 x 16 lattice, for which we do not consider 
the full line shape obtained with the Max-Ent method 
to be stable due to larger statistical errors than for the 
smaller systems. The first two cumulants can neverthe- 
less be estimated. The Max-Ent and series results agree 
very well for the larger systems. The exact results for 
the non-square lattices do not show a regular size depen- 
dence, whereas the L x L lattices do. The moments for 
the L x L lattices increase monotonically with L. How- 
ever, there is a clear maximum in the second cumulant 
for L = 6. This is likely caused by the lack of weight be- 
tween the two dominant peaks for this lattice size. With 
growing size the gap should gradually be filled in by other 
peaks, leading to a decreasing second cumulant. Judging 
from Fig. ||, the results for the largest systems (16 x 16 
for Mi and M-x and 10 x 10 for M3) should represent 
the thermodynamic limit within statistical errors. We 
then have M x = 3.59 ± 0.01, M 2 = 0.79 ± 0.03, and 
M 3 = 0.95 ± 0.08. 

We now return to the line shape. The Max-Ent spec- 
tra displayed in Fig. |^ show a considerable dependence 
on the lattice size. The trend for L > 6 appears to be the 
development of a well defined main peak at uj/J ~ 3.5, 
as well as some strengthening of the tail up to uj/J ~ 7. 
Comparing with the spinwave results for the two-magnon 
profiles shown in Fig. 0, the 10 x 10 Max-Ent spectrum 
is clearly much broade^than the narrow peak obtained 
by Canali and Girvin,tZl but not quite as much broad- 
ened towards lower frequencies as the Chubukov-Frenkel 
profilec.3 obtained by setting S = 1/2 at a later stage 
of the calculation. Since the Max-Ent method can be 
expected to cause some broadening and the trend with 
increasing the lattice size appears to be a narrowing of 
the dominant peak, we conclude that the actual peak in 
the thermodynamic limit should be narrower than that 
obtained by Chubukov and Frenkel. 

In the exact 6x6 result there are contributions in the 
frequency range u>/ J ss 4.5 — 7 which are not present in 
the spinwave result for the same lattice (see Fig. 0). This 



weight is therefore most likely dominated by processes 
involving more than two magnons. The Max-Ent result 
for the 10 x 10 lattice also shows a tail extending up to 
lu/J f» 7. The total weight above the spinwave theory 
two-magnon cut-off lu/J = 4.63 does, however, remain 
at w 10 — 15%, as previously argued on the basis of-the 
1/S 2 spinwave results and the frequency cumulants. EH 

Canali and Girvin argued that the two-magnon pro- 
file is very little affected by the higher-order processes, 
and that the four-magnon contribution should-bc a peak 
well separated from the two-magnon profile.EZI We now 
consider an approach to testing this hypothesis numer- 
ically, independently of the Max-Ent method. We as- 
sume a spectrum consisting of the Canali-Girvin two- 
magnon profile P(lu) shown in Fig. ^, and a Gaussian 
G ai (uj — L04) of width (74 centered at to = 104 for modeling 
the higher-order contribution. In order to account for 
a possible further frequency shift, we use a phcnomeno- 
logical frequency renormalization Z in the two-magnon 
profile. The full spectrum is hence 



I{uo) = A 2 P(Zlu) + A 4 G a 



(oj- 



- w 4 ), 



(31) 



where P and G CT4 are both normalized to one, and hence 
A2 + A4 = 1. We then have four parameters; Z, A4, 
uj 4. , and o~4 which can be adjusted to give the best con- 
sistency with the imaginary-time data. Note that P(lu) 
already contains the spinwave renormalization factor to 
order 1/5 2 , and hence our Z should be close to 1 for this 
treatment to be consistent. 

For the 10 x 10 lattice, the imaginary-time data can 
indeed be very well accounted for by this spectrum. We 
obtain the parameters Z = 0.97, A4 = 0.40, 104 = 4.1 and 
0-4 = 1.1. The resulting spectrum is shown in Fig. ^. The 
data for 16 x 16 spins can also be very well fit to the form 
considered, and the parameters are not changed much 
from the 10 x 10 ones. This spectrum is also shown in 
Fig. [| The parameters of the Gaussian are such that it 
is completely merged with the two-magnon profile. This 
is clearly consistent with both the 6x6 Lanczos and the 
Max-Ent results, which did not show any significant gap 
between the main peak and the high-frequency weight. 
In Fig. H the weight of the Gaussian also extends to the 
low-frequency side of the two-magnon peak, and therefore 
has the effect of broadening it. Therefore, the relative 
weight of about 40% of the secondary peak cannot be in- 
terpreted directly as the total four-magnon (and higher) 
contribution, but also likely reflects that the two-magnon 
profile from spinwave theory is too narrow. We have also 
carried out fits to two Gaussians, and then find that the 
dominant one is at a position « 3.2J, and the second 
one again is at ~ 4 — 4.5 J. However, the uncertainty in 
the width of the dominant peak is large, and therefore 
this method cannot be used to accurately determine the 
width. Based on the other approaches we have discussed, 
we can nevertheless conclude that the standard spinwave 
two-magnon profile is too narrow, but by how much is 
not completely clear. The profile shown in Fig. || likely 
represents a lower bound of the width. 
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We also attempted a similar fitting procedure using the 
Chubukov-Frenkel two-magnon result as the dominant 
feature. However, we found that it was not possible to 
obtain any good fit to the QMC data in this case, due to 
the, apparently, too high low-energy weight. 

In Fig. ^ we also show an experimental spectrum for 
La2CuC>4, with the frequency scale adjusted to give the 
same peak position ui/ J w 3.25 as the QMC-spinwave fit. 
This peak position corresponds to an exchange constant 
J = 1440 K for the experimental system, which is in good 
agreement with J ss 1500 K obtained from Neutron scat- 
tering and NMR experiments. Although the experimen- 
tal spectrum is somewhat broader than our result, there 
is a quite good agreement with the distribution of the 
weight present above the two-magnon cut-off frequency. 
Since the width of the theoretical spectrum shown here 
most probably is a lower bound of the actual width, we do 
not consider the deviations from the experimental spec- 
trum serious. Comparing with the two-magnon spinwave 
spectra shown in Fig. |^, it is clear that our present fit- 
ted spectrum is considerably closer to the experimental 
result. As will be discussed further in Sec. VI, the width 
of the peak is such that the further broadening required 
to match the experimental spectrum could quite easily 
be achieved by spin-phonon couplings, as has been sug- 
gested by several groups. 

V. FINITE-TEMPERATURE RESULTS 

In this section we present results of QMC and Max- 
Ent calculation^ carried out at temperatures T j J = 0.25, 
0.5, and 1.0. c3 Raman spectra for a 4 x 4 lattice at 
these temperatures were previously obtained by Bacci 
and Gagliano using exact diagonalization.Ej Recently, fi- 
nite temperature Lanczos calculations for lattices with up 
to 20 sites were presented by Prelovsek and JaklicO Here 
we compare QMC+Max-Ent results for systems with 4x4 
and 16 x 16 spins. The latter size should be sufficient for 
obtaining thermodynamic limit results at the tempera- 
tures considered. 

Fig. ^ shows the imaginary-time correlation functions. 
For the temperatures considered here, g(r) can be accu- 
rately evaluated for the whole range < t < (3. The 
slower decay with r for the larger lattice indicates the 
presence of more low-frequency weight as the system size 
increases. This is confirmed by the Max-Ent results for 
the real-frequency spectra, shown in Fig. [IT]. The results 
for 4 x 4 spins are in reasonable agreement with exact 
diagonalization results if one includes some rather large 
broadening of the 5- functions. In Fig. [ll] we have graphed 
the exact results as histograms, with the bin width for 
each temperature chosen large enough to remove most, 
but not all, of the jagged structure due to the discrete 
finite-size spectrum. It is clear that the Max-Ent method 
cannot capture the fine-structure of the spectrum, and 
instead gives a single rounded shape. Nevertheless, the 



region of dominant spectral weight and its temperature 
variations are well reproduced. The low-frequency peak 
in the exact 4x4 spectra at high temperatures is due 
to degeneracies present for this small latticeEj (i.e., the 
peak is actually at uj = 0). 

Our 16 x 16 results show a faster enhancement of the 
low-frequency spectral weight as the temperature is in- 
creased above T/J w 0.25. This difference between lat- 
tice sizes is likely due to the presence of large finite-size 
gaps in the level spectrum of the 4x4 system. Naturally, 
as T — > oo the system size dependence should diminish, 
and this is seen already at T — 1.0 in Fig. [ll]. The finite- 
temperatupe, spectra calculated for 20 sites by Prelovsek 
and Jaklicca for T/J = 1.0 and 0.5 are in reasonable 
agreement with our 16 x 16 results, again taking into ac- 
count a Max-Ent broadening of our spectra. However, at 
T/J = 0.5, judging from the rather large differences be- 
tween the exact results for N = 16 (Ref. |o|) and N = 20 
(Rcf. and the slow approach to the thermodynamic 
limit discussed in sec. IV, it is likely that the N = 20 
spectrum has not yet converged to its infinite-size shape. 
The actual width at this temperature should therefore 
be something intermediate between our 16 x 16 Max-Ent 
result and the previously obtained N = 20 profile. 

Experimentally, spectra taken at room temperature do 
not differ sig nific antly from ones obtained at very low 
temperatures.li3liil As the temperature is elevated to to 
T/J ps 0.5 thete is a significant increase in the weight 
below oj 2 J.EJ This feature is indeed quite well repro- 
duced by our result for 16 x 16 spins. 

The spectra shown in Fig. [ll] are all normalized to 1 . 
The temperature dependence of the integrated intensity 
is a quantity of experimental interest. We define two 
intensities: 

/oo 
dujA(uj), (32a) 
-OO 

I 2 = duA{uj). (32b) 
Jo 

These definitions are equivalent at T = 0, but differ at 
finite T due to spectral weight at negative frequencies, 
with A{— w) = e~@ UJ A(uj). I\ can be directly obtained 
from the imaginary-time data as G(t = 0), whereas I 2 
is calculated by integrating the real frequency spectrum 
obtained using the ME method. Fig. |l2| shows both in- 
tensities vs. T for 4x4 and 16 x 16 lattices. Up to 
T/J w 0.25, 1\ « I2, owing to the absence of significant 
low-frequency weight at these temperatures. At higher 
temperatures l\ > I2, but even at T/J sa 0.5 the dif- 
ference is small. For the 4x4 system the intensity I2 
increases by ~ 14% as the temperature is decreased from 
T/J = 0.5 to T/J = 0, and for 16 x 16 by « 9%. 
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VI. SUMMARY AND DISCUSSION 

We have presented numerical results for the Bi g spec- 
trum of the Heisenberg model within Loudon-Fleury the- 
ory. We obtained Lanczos exact diagonalization results 
for up to 6 x 6 spins, and carried out QMC simulations 
for up to 16 x 16 spins. We compared the results with 
spin wave theory. Our main results and conclusions are 
the following: 

1) Comparing spinwave theory and exact diagonaliza- 
tion results for the same lattice sizes, we find that for a 
given cluster the number of dominant peaks is the same 
in both cases. However, both the positions of the peaks 
and their relative weights are different. Most notably, 
for the 6x6 lattice there are two dominant peaks, the 
separation of which is 1.5 times larger in the exact result. 
Assuming that the trend persists for larger lattices, this 
indicates that spinwave theory underestimates the width 
of the dominant B\ g peak in the thermodynamic limit. 

2) Our results of Maximum-Entropy analytic contin- 
uation of QMC imaginary-time data is also consistent 
with a peak width larger than that of the spinwave two- 
magnon peak. The first three frequency cumulants are in 
excellent agreement with previous results of a series ex- 
pansion around the Ising model. We estimate the cumu- 
lants in the thermodynamic limit to be M\ = 3.59±0.01, 
M 2 = 0.79 ± 0.03, and M 3 = 0.95 ± 0.08. 

3) In order to test the 1/S 2 spinwave theory predic- 
tion of a four-magnon jpppfile well separated from the 
main two-magnon peakJlH we carried out a fit of the 
QMC imaginary-time data to a spectrum consisting of 
the spinwave two-magnon peak and a Gaussian at higher 
frequency. We found that this type of spectrum indeed 
describes the data well. The fitted Gaussian is cen- 
tered at ui/J sa 4.1, and is so broad that it is com- 
pletely merged together with the two-magnon structure 
peaked at uj/J w 3.25. The resulting spectrum resem- 
bles a typical experimental B\ g profile for La2Cu03 with 
an exchange J k, 1400 K. The experimental peak is still 
slightly broader, but there is a considerable improvement 
in comparison with the standard spinwave theory two- 
magnon profile. 

4) The imaginary-time data cannot be fitted using the. 
two-magnon profile obtained by Chubukov and Frenkelcfl 
by expanding their spinwave spectrum around its peak 
position before setting S = 1/2 in the calculation. This 
is due to the significantly stronger low-frequency weight 
present in this spectrum. 

5) At finite temperature we find a significant increase 
in spectral weight below oj w 2 J for T/J > 0.25, in 
agreement witb-|experimental results for antiferromag- 
netic cuprates.EJ We also find that this effect is sup- 
pressed in the 4x4 system, due to the finite-size gaps. 
The temperature dependence of the integrated scattering 
intensity is weak. 

Our results hence confirm that LF theory can account 
for some of the main features of typical B\ g spectra 



observed experimentally for antiferromagnetic cuprates 
such as La2Cu04. Our new evidence for a profile sig- 
nificantly broader than that obtained in spinwave theory 
support in part the early claim by Singh et alSn that the 
broadening is due to the strong quantum fluctuations of 
the Heisenberg model with 5 = 1/2 (note that spinwave 
theory is in good agreement, with experimental results 
for quasi-2D S = 1 systemscj). However, typical exper- 
imental spectra are still broader, and extend to slightly 
higher frequencies. The shoulder-like feature observed in 
some experiments at lu 4J is also not present in our 
results, although we find evidence that the four-magnon 
contribution has its maximum in this regime. Hence, 
although our results show a better agreement with ex- 
periments than, prev ious numerical results obtained for 
smaller lattices JIM2I the Heisenberg-LF mechanism does 
not appear to fully account for the experimental Raman 
scattering, as has been noted in several previous studies. 
The fact that there is no A\ g scattering within this theory 
of course also implies that other additional mechanisms 
have to be active. 

Chubukov and Frenkel recently suggested that reso- 
nant processes not contained within LF theory are im- 
portant in typical experiments, for which the frequency 
cjj n of the incident light is coraflarable to the charge trans- 
fer gap of the Cu02 planes .Eil We agree that resonance 
effects are most likely needed to explain the dependence 
of the total scattering intensity and the line shape on 
Wi n jL3 but note that the dominant features of the profile 
do not show much-dependence on uj- ln for most of the fre- 
quencies studied S3l3 Based on the improved agreement 
with experiments obtained here within LF theory, we be- 
lieve that the main features of the B\ g spectrum are due 
to the LF mechanism, and that a further broadening of 
the spectrum could be achieved by magnon damping due 
to phonons. 

Motivated by experiments carried out at high temper- 
atures, Knoll et al. suggested that spin-lattice interac- 
tions may be r.esponsible for the broadening of the Ra- 
man spectrum.uJ Spinwave calculations including a phe- 
nomenologipal magnon life-time give some support to 
these ideaslL!l Several different calculations explicitly in- 
cluding jH agno n-phonon coupling have been presented 
recentlyE3Scia Using an adiabatic approximation for the 
phonons leads to a Heisenberg model with random cou- 
pling constants. Numerically studying such random lat- 
tices with 4x4 spins and assuming a standard LF cou- 
pling, Nori et al. found that the B\ g spectrum can indeed 
be broadened by this mechanism, and that also A\ g scat- 
tering can become significant H3 However, in this calcu- 
lation, the strength of the randomness required in order 
to reproduce the width of the experimental B\ g spec- 
trum appears to be rather large (using a Gaussian dis- 
tribution for the nearest-neighbor couplings Jij , a width 
a w 0.5(Jij) was required)M3 Nori et al. argued that 
such strong disorder can be caused by incoherent atomic 
displacements. Nevertheless, in the absence of other evi- 
dence for the presence of large fluctuations in the Heisen- 
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berg couplings, it would be desirable to reproduce the 
broadened spectrum with a narrower coupling distribu- 
tion. 

One reason for the strong disorder required in the cal- 
culation of Nori et al. could be the small size of the lattice 
used.E3 As we have seen, the pure 4x4 system only has 
a single dominant two-magnon 5-function at u> —JL98J, 
and two weaker peaks at to w 4.5 J and lo w 5.5 J S3 It is 
clear that a considerably weaker disorder would suffice to 
broaden the spectrum if one starts from the much broader 
pure-system profile obtained here for larger lattices. 

The type of QMC and Max-Ent calculations presented 
here could in principle be carried out also for disor- 
dered spin, systems, and even including fully dynamic 
phonons.EJ The suggested effects of magnon-phonon cou- 
pling could hence be investigated more rigorously than 
previously, using larger lattices. Although a^ecent exact 
diagonalization study by Reilly and RojoE3 give some 
support for the validity of an adiabatic approximation 
for the phonons, calculations with full dynamic phonons 
should also be carried out for larger lattices. Limits 
on the strength of the phonon-magnon coupling (or the 
width of the disorder distribution in the adiabatic ap- 
proach) could be established by carrying out QMC cal- 
culations of, e.g., the temperature dependence pfjthe spin 
correlation lengthB and NMR relaxation ratesE3 for sys- 
tems including lattice vibrations or static disorder. 
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APPENDIX A: QMC CALCULATIONS OF THE 
RAMAN CORRELATION FUNCTIONS 

In this Appendix we describe the calculation of 
the imaginary-time correlation function (KS^with the 
stochastic series expansion (SSE) methodElEj In order 
to reduce the statistical fluctuations, we use the spin- 
rotational invariance of the Heisenberg Hamiltonian to 
construct an estimator less noisy than the obvious one. 
We also derive direct estimators for the r-derivatives of 
G(r). In order to establish the notation, we first very 
briefly outline the formalism of the SSE algorithm. More 
details of the implementation of thifpaon-standard gener- 
alization of Handscomb's methocSo for the 2D Heisen- 
berg model can be found in Ref. 09. 



In order to apply the SSE technique, the Hamiltonian 
is first written as 



2N 



ft - ~^ ]T) - ft 2 ^ + ~0~> 



(Al) 



6=1 



where 6 is a link connecting a pair of nearest-neighbor 
sites (i(b),j(b)), and the operators Hij, and H 2 £ are de- 
fined as 



Ai,b = 2[| 



S i(b) S j(b)\ 



(A2a) 
(A2b) 

An exact expression for an operator expectation value 



H 2 ,b - St(b) S j(b) 



S i(b) S j(b)- 



(A) = -Tr{ie-^}, 



Tr{e"^} 



(A3) 



at inverse temperature (3 — J/T, is obtained by Tay- 
lor expanding exp(—(3H) and writing the traces as sums 
over diagonal matrix elements in the basis {|a)} = 
{\Slj-j . , Sff)}. The partition function then takes the 
formEJ 



-EEE^RHII^, 

a n S„ 1 = 1 



(A4) 



where S n is a sequence of index pairs defining the oper- 
ator string nr=i H aiM-> 



S n = [ai,b 1 ][a 2 ,b 2 ] ■ ■ ■ [a n ,b n ], 



(A5) 



with <n £ {1)2}, bi € {1, ...,2iV}, and n 2 denotes the 
total number of index pairs (operators) [a.j, bi] with a, = 
2. Both H\ y b and H 2> b can act only on states where the 
spins at sites i(b) and j(b) are antiparallel. leaves 
such a state unchanged, whereas H 2 ^ flips the spin pair. 
Defining a propagated state 



\a{p))=i[H, 



<*{,&! I 



*(0)> = 



(A6) 



i=i 



a contributing (a,S n ) must clearly satisfy the periodic- 
ity condition |a(n)) = |a(0)). In an allowed sequence 
S n , the links b corresponding to the spin- flipping opera- 
tors [2, b] present must therefore form only closed loops. 
For a lattice with L x L sites and L even, this implies 
that the num ber n 2 must be even, and hence that all 
terms in Eq. (A4) are positive and can be used as rel- 
ative probabilities in a Monte Carlo algorithm (this is 
true for any non- frustrated system). Since any non-zero 
matrix element in (A4) is equal to one, the weight factor 
corresponding to a contributing (a, S n ) is simply given 
by 



W(a,S n ) 



(0/2) 



(A7) 
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The algorithm for sampling the configurations (a, S n ) is 
described in Ref. |9| 

In order to obtain an expression for G(r) in terms of 
the states |a(p)) and the index sequence S n used in the 
simulation, the expectation value is first written in terms 
of the operators H a ,b as 



G(r) 



E E 

ai,a 2 b±,b 2 



g:±(t), 



a 2 ,b. 



where 



G a a ±{r) 



(H a2M ( T )H aiM (0)), 



(A8) 



(A9) 



and B\g symmetry corresponds to Pb lt b 2 — 1 f° r links b\ 
and b 2 which are parallel to each other, and Pb x ,b 2 = — 1 
for perpendicular links. P roce eding as in the derivation 
of the partition function (A.4), the exponentials in the 
expression 



a 2 ,b 

z ^ 



(A10) 



are Taylor expanded and all powers of H are written 
as sums of products of the operators H a ,b- There is 
then a one-to-one correspondence between the terms in 
(r) and Eq. (|A4|). Dividing out the factor corre- 



Tf<ll,&l 

T a 2 ,b^ 



sponding to the configuration weight (A7) gives the av- 
erage in the form of a function of S n : 



' n-2 



F(r,n-,m)N^jHrn) 



\m=0 



where 



F(r,', 



and N* 1 ' 6 * (m) is the number of times the operators 
[ox, bi] and [d 2 , b 2 ] occur in S n (in the given order) sepa- 
rated by m other operators. Hence, measuring G^'^r) 
simply amounts to finding all pairs of operators [ax,&x] 
and \a 2 , 6 2] in the sequence S n . The contribution to 
Eq. ( A.10 ) of each pair is a function o f the relative sepa- 
ration of the operators, given by Eq. (A.12). 

In order to obtain a simple expression for the full cor- 
relation function G(t) it is useful to introduce a function 
X(p), such that X(p) = +1 if the p:th operator in S n acts 
on a link in the cc-direction, and X(p) = —1 if it acts on 
a y link. Numbering the bonds such that < b < N 
correspond x-bonds, and N + 1 < b < 2N correspond to 
y-bonds, the definition is hence 



! (n- 1)! 



m 



2)!m! 



(All) 



(A12) 



X(p) = 



b p <N 
b p >N. 



(A13) 



Eqs. (|A9|) and (All) then give 



G ( r ) = ( E E F ( r > m ~ ^ x (p) x (p + ■ 



(A14) 

where of course X(p) is perio dic; X (n + 1) = X(l). 

In practice the estimator ( |A14| ) is rather noisy. An 
improved estimator can be constructed as follows. First, 
the function X{p) is written as a sum of two terms, 



X(p)=X 1 ( P )+X 2 (p), 



(A15) 



where X t (p) = ±1 (t — 1, 2) for x and y bonds, as before, 
if the p:th operator in S n , [a p , b p ] , has a p = t, but X t (ja) = 
if a p ^ t. Hence 



Xt( P ) 



, a p ^ t 



b p <N 
b p >N 



(A16) 



If a p = 1, Xi(p) can be averaged over all 2N choices 
of operators [1,6] at position p. The weight W(a, S' n ) 
corresponding to a sequence S' n obtained by replacing the 
current operator [1, b p ) at p in S n is equal to the current 
weight W(a, S n ) if the corresponding spins at sites i(b) 
and j(b) are antiparallcl in the propagated state \a(p)), 
and is zero otherwise. Hence, Xi(p) can be redefined as 

X 1{ P) = \ [N * {P) - N J {P)] I 2N (AIT) 

where N^(p) is the number of antiparallel nearest- 
neighbor spin pairs in the 7-direction in \a(p)). One can 
easily verify that this averaged estimator can be used in 
products with both X\ and X 2 . H ence, improved es- 
timators for X(p)X(p + m) in Eq. ( A14 ) can be used 
for the terms X\X\, X1X2, and X 2 Xi. For X 2 (p) no 
simple re-definition in terms of single-operator averaging 
can be constructed (replacing a single operator [2, b] with 
any other operator always leads to a non-contributing 
term), and hence the X 2 X 2 contribution to (A14) re- 
mains noisy. However, the rotational invariance of the 
Heisenberg Hamiltonian implies that, 

(X 2 (p)X 2 (p + m)) = (2X 1 {p)X 1 (p + m)) + 

\(X 1 {p)X 2 {p + m) + X 2 (jp)Xi(p + to)), (A18) 

and therefore the X 2 X 2 term does not even have to be 
evaluated. The final result for the improved estimator 
for G(t) is hence 

/ ^ n n— 1 



p—1 m=l 



Xi {p)X 2 (p + m)+ X 2 (p)X 1 (p + m) 



(A19) 



It should be noted that the function F{j, n; m) is sharply 
peaked around m ps jit / (3 for large f3, so that typically 
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only a small fraction of the terms in (A19) actually have 
to be e valua ted. 

Eq. ( |A19|) is valid for any < r < /3, and the r- 
dependence appears only in the function F(r,n;m). In 
contrast to standard Trotter-based QMC methods, the 
method discussed here can therefore be used to directly 
calculate also r-derivatives of imaginary-time dependent 
correlation functions. An expression for the n:th deriva- 
tive of G(r) is simply obtained by replacing F in (A19) 
by its n:th derivative: 



4 P. Bourg es, H. Casalta, A. S. Ivanov, and D. Petitgrand, 
preprint, |sond-mat/9708060| (1997). 



G (n) (i 



d n G(r) 
dr n 



7i n—1 
p— 1 m— 1 



'd n F(T, n, m — I) 



X 1 (p)X 2 (p + m) + X 2 (p)X x (p + m) 



(A20) 



As discussed in Sec. IV, derivatives can be used as supple- 
mentary information in a numerical analytic continuation 
to real frequency. The derivatives at t = are of special 
interest, as they are related to moments of the spectral 
function [see Eq. (|29|)]. 

We end this Appendix with a demonstration that the 
simulation results for G(t) are indeed free from system- 
atic errors. Since the absolute Raman scattering intensity 
is not contained in the LF theory, the amplitude of /(w), 
and hence of G(r), is irrelevant, and instead of G(r) one 
can consider the ratio 



g{r) = G(t)/G(0). 



(A21) 



Fig. 13 shows the QMC result for this quantity calculated 
on a 4 x 4 lattice, along with the exact result obtained 
from I(oj) calculated using exact diagonalization. The 
statistical error of the QMC result is in the fifth decimal 
digit, and there is excellent agreement with the exact 
result within this accuracy. The absence of detectable 
systematical errors in the QMC result for g(r) is hence 
confirmed. Since g(r) decays exponentially, the relative 
statistical error grows rapidly with r, and for r > 3 accu- 
rate results can not be easily obtained. This is the case 
also for larger systems. 
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FIG. 1. Spinwave theory results for the Bi g Raman profile 
calculated on small lattices with L x L sites. The dashed 
curves are the results with the interactions neglected, and the 
solid ones are with interactions in the final states included 
at the RPA level. A damping e = 0.05J has been used to 
broaden the (^-functions. 



FIG. 2. Spinwave theory results for the B\ g two-magnon 
profile in the thermodynamic limit, compared with the ex- 
perimental spectrum for La2Cu04 discussed in Ref. UfA (bold 
solid curve). The solid curve corresponds to Eq. ( |18| ) with 
a spinwave renormalization factor Z c — 1.158. The dotted 
curve is the result by Canali and Girvin, which includes also 
quantum fluctuations in the ground state. The dashed line 
is the result by Chubukov and Frenkel, obtained by further 
expanding the line shape (|l|) in 1/S before setting 5 = 1/2. 
All curves are normalized to one. The frequency scale of the 
experimental spectrum has been adjusted to give a peak po- 
sition in rough agreement with the theoretical curves. 



FIG. 3. Exact diagonalization results for the Bi g spectrum 
for different small lattices with N sites (solid curves). The 
dashed curves are the corresponding RPA-spinwave results. 
The 5-functions of the exact results have been broadened us- 
ing a damping e = 0.1J, and all the spectra ar normalized to 
one. The spinwave results have been given a smaller damping 
and a different normalization in order to more clearly show 
the peak positions. 



FIG. 4. QMC results for In [g(r)] (upper panel) for different 
system sizes, and the relative statistical errors of <?(r) (lower 
panel) . 



FIG. 5. Results of Max-Ent analytic continuation of 10 
bootstrap samples of QMC imaginary-time data generated 
for a 4 x 4 system. 



FIG. 6. Bootstrap-averaged Max-Ent results for the B\ g 
spectrum for different lattices (solid curves). The 4x4 and 
6x6 results are compared with the corresponding exact diag- 
onalization results with a damping e/J = 0.1 (dotted curves). 



FIG. 7. QMC results for the short-time behavior of G (1) /G 
(upper panel) and G' 2 ' / G (lower panel) for systems of linear 
sizes L = 4 (solid circle), 6 (open circles), 8 (solid squares), 
and 10 (open squares). The solid curves are obtained from 
the Max-Ent analytic continuation. 
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FIG. 8. The first three frequency cumulants of the 
Max-Ent spectra vs. the inverse system size (solid circles with 
error bars). The open circles are the exact diagonalization 
results. The previous infinite-size riesults from a series ex- 
pansion, calculated by Singh et al.f^ are indicated by the 
horizontal dashed lines (result ± estimated error). 

FIG. 9. Big spectrum obtained by a fit of imaginawi-time 
QMC data to the Canali-Girvin two-magnon profildij plus 
a Gaussian. The almost indistinguishable solid and dashed 
curves are for a 10 x 10 and a 16 x 16 lattice, respectively. 
The bold curve is the experimental spectrum for La2Cu04 
discussed in Ref. with the frequency scale adjusted to give 
the same peak position as the theoretical results (correspond- 
ing to an exchange J = 1440 K). 

FIG. 10. The logarithm of the normalized imaginary-time 
correlator g(r) vs. r for 4x4 (dashed curves) and 16 x 16 
(solid curves) lattices at different temperatures. 

FIG. 11. Max-Ent results for the Bi g spectrum of 4 x 4 
(dashed curves) and 16 x 16 (solid curves) lattices at different 
temperatures. The histograms represent the exact results for 
the 4x4 lattice. 

FIG. 12. Integrated B\ g scattering intensities vs. temper- 
ature for 4x4 (open symbols) and 10 x 10 (solid symbols) 
lattices. Circles are for Ii (using all frequencies), and squares 
for I2 (using positive frequencies only). 

FIG. 13. Upper panel: QMC results for g(r) = G(r)/G{0) 
of a 4 x 4 system at inverse temperature /3 = 32 (solid circles), 
compared with the exact ground state result (solid curve). 
The inset shows the regime 1.75 < r < 3 on a more detailed 
scale. Lower panel: The deviation of the QMC data from the 
exact result, multiplied by 10 4 . The dashed curves indicate 
the statistical errors. 
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Fig. 2, A. W. Sandvik et al. 




Fig. 3, A. W. Sandviket al. 




Fig. 4, A. W. Sandvik et al. 



15 




Fig. 5, A. W. Sandvik et al. 
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